#### SCRIPT FOR "How Voters Respond to Presidential Assaults on Checks and Balances: Evidence from a Survey Experiment in Turkey"
#### Comparative Political Studies, 2021


#### SETUP ####
setwd("~/Dropbox/Winning prospects and changing preferences/Publication") 
# Note: Change the working directory to the folder to which you downloaded the data. 
rm(list = ls()) 

## Packages
need <- c("foreign", "stargazer", "arm", "openxlsx", "clubSandwich", "nnet", "margins", "readstata13", "sjPlot", 
          "sjmisc", "ggplot2", "lmtest", "multiwayvcov", "DescTools") 
have <- need %in% rownames(installed.packages()) 
if(any(!have)) install.packages(need[!have]) 
invisible(lapply(need, library, character.only=T)) 
rm(have, need)


## Opening the datasets
load("ExecAggrCitizenResponse_Data.Rdata")

## Dividing the datasets
execaggr_yes <- execaggr[(execaggr$opposition_narrow == 0 & is.na(execaggr$opposition_narrow) == 0),]
execaggr_no <- execaggr[(execaggr$opposition_narrow == 1 & is.na(execaggr$opposition_narrow) == 0),]

#### DESCRIPTIVE STATISTICS ####
# Descriptive statistics for outcome variables [Table 1]
stargazer(execaggr_yes[c("decree_power", "budget_power", "system_endorsement", "executive_president", "noreturn_reduced")], 
          type="text")
stargazer(execaggr_no[c("decree_power", "budget_power", "system_endorsement", "executive_president", "noreturn_reduced")], 
          type="text")

t.test(execaggr_yes$decree_power[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$decree_power[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$budget_power[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$budget_power[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$system_endorsement[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$system_endorsement[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$executive_president[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$executive_president[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$noreturn_reduced[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$noreturn_reduced[execaggr_yes$treatment_erdoganloss==1])

t.test(execaggr_no$decree_power[execaggr_no$treatment_erdoganloss==0], execaggr_no$decree_power[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$budget_power[execaggr_no$treatment_erdoganloss==0], execaggr_no$budget_power[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$system_endorsement[execaggr_no$treatment_erdoganloss==0], execaggr_no$system_endorsement[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$executive_president[execaggr_no$treatment_erdoganloss==0], execaggr_no$executive_president[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$noreturn_reduced[execaggr_no$treatment_erdoganloss==0], execaggr_no$noreturn_reduced[execaggr_no$treatment_erdoganloss==1])

# Descriptives of the moderator variables (Table 3 in the paper)
stargazer(execaggr_yes[c("anxiety_economy", "soc_polar1", "soc_polar2", "soc_polar3", "socpol_additive")], 
          type="text")
stargazer(execaggr_no[c("anxiety_economy", "soc_polar1", "soc_polar2", "soc_polar3", "socpol_additive")], 
          type="text")

t.test(execaggr_yes$anxiety_economy[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$anxiety_economy[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$soc_polar1[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$soc_polar1[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$soc_polar2[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$soc_polar2[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$soc_polar3[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$soc_polar3[execaggr_yes$treatment_erdoganloss==1])
t.test(execaggr_yes$socpol_additive[execaggr_yes$treatment_erdoganloss==0], execaggr_yes$socpol_additive[execaggr_yes$treatment_erdoganloss==1])

t.test(execaggr_no$anxiety_economy[execaggr_no$treatment_erdoganloss==0], execaggr_no$anxiety_economy[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$soc_polar1[execaggr_no$treatment_erdoganloss==0], execaggr_no$soc_polar1[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$soc_polar2[execaggr_no$treatment_erdoganloss==0], execaggr_no$soc_polar2[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$soc_polar3[execaggr_no$treatment_erdoganloss==0], execaggr_no$soc_polar3[execaggr_no$treatment_erdoganloss==1])
t.test(execaggr_no$socpol_additive[execaggr_no$treatment_erdoganloss==0], execaggr_no$socpol_additive[execaggr_no$treatment_erdoganloss==1])


#### ANALYSES HYPOTHESIS 1 ####
## Creating the factor variable for the treatment first
execaggr$tretman <- ifelse(execaggr$treatment_erdoganloss==1, "2 Incumbent Loses", "1 Incumbent Wins")
execaggr_yes$tretman <- ifelse(execaggr_yes$treatment_erdoganloss==1, "2 Incumbent Loses", "1 Incumbent Wins")
execaggr_no$tretman <- ifelse(execaggr_no$treatment_erdoganloss==1, "2 Incumbent Loses", "1 Incumbent Wins")

# Models
fit.h1.decree.yea <- glm(decree_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h1.decree.nay <- glm(decree_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))
fit.h1.budget.yea <- glm(budget_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h1.budget.nay <- glm(budget_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))
fit.h1.system.yea <- polr(factor(system_endorsement) ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_yes, Hess=T)
fit.h1.system.nay <- polr(factor(system_endorsement) ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_no, Hess=T)
fit.h1.partis.yea <- glm(executive_president ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h1.partis.nay <- glm(executive_president ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))
fit.h1.noretu.yea <- glm(noreturn_reduced ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h1.noretu.nay <- glm(noreturn_reduced ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

# Clustering-robust standard errors
se.h1.decree.yea <- list(coef_test(fit.h1.decree.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h1.decree.nay <- list(coef_test(fit.h1.decree.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.h1.budget.yea <- list(coef_test(fit.h1.budget.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h1.budget.nay <- list(coef_test(fit.h1.budget.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.h1.partis.yea <- list(coef_test(fit.h1.partis.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h1.partis.nay <- list(coef_test(fit.h1.partis.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.h1.noretu.yea <- list(coef_test(fit.h1.noretu.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h1.noretu.nay <- list(coef_test(fit.h1.noretu.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])

# Block-bootstrapped standard errors for the ordinal logistic regression: Column 5
est <- function(data1){
  fit.h1.system.yea <- polr(factor(system_endorsement) ~ tretman + rte_firstround + knowledge + left_right, data=data1, Hess=T)
  results <-  c(fit.h1.system.yea$coefficients)
  return(results)
}

resamp.block <- function(data){
  lookup <- split(1:nrow(data), data$ilcemah)
  mahalleler <- names(lookup)
  mahalleler.drawn <- sample(mahalleler, length(mahalleler), replace = TRUE)
  star <- unlist(lookup[mahalleler.drawn])
  samp.blo.data <- data[star,]
}

ests <- matrix(NA, 4, 5000)
set.seed(34718)
for(i in 1:5000){
  tryCatch({
    data2 <- resamp.block(execaggr_yes)
    ests[,i] <- est(data2)
  }, error=function(e){})
  if(i %% 10 == 0) cat(i," ")
}

mns <- c(mean(ests[1,], na.rm=T), mean(ests[2,], na.rm=T), mean(ests[3,], na.rm=T), mean(ests[4,], na.rm=T))
sds <- c(sd(ests[1,], na.rm=T), sd(ests[2,], na.rm=T), sd(ests[3,], na.rm=T), sd(ests[4,], na.rm=T))
tscs <- mns/sds 
mns
sds <- list(c(sds[1], sds[2], sds[3], sds[4])) # These are used as the standard errors in Columns 5. 
sds
tscs
(1-pt(abs(tscs), 500))*2

# Block-bootstrapped standard errors for the ordinal logistic regression: Column 6
est2 <- function(data1){
  fit.h1.system.nay <- polr(factor(system_endorsement) ~ tretman + rte_firstround + knowledge + left_right, data=data1, Hess=T)
  results <-  c(fit.h1.system.nay$coefficients)
  return(results)
}

ests2 <- matrix(NA, 4, 5000)
set.seed(34719)
for(i in 1:5000){
  tryCatch({
    data2 <- resamp.block(execaggr_no)
    ests2[,i] <- est2(data2)
  }, error=function(e){})
  if(i %% 10 == 0) cat(i," ")
}

mns2 <- c(mean(ests2[1,], na.rm=T), mean(ests2[2,], na.rm=T), mean(ests2[3,], na.rm=T), mean(ests2[4,], na.rm=T))
sds2 <- c(sd(ests2[1,], na.rm=T), sd(ests2[2,], na.rm=T), sd(ests2[3,], na.rm=T), sd(ests2[4,], na.rm=T))
tscs2 <- mns2/sds2 
mns2
sds2 <- list(c(sds2[1], sds2[2], sds2[3], sds2[4])) # These are used as the standard errors in Columns 6. 
sds2
tscs2
(1-pt(abs(tscs2), 500))*2

pseudo.r2s.h1 <- c("Psuedo R-sq.", 
                round(PseudoR2(fit.h1.decree.yea, "McFadden"), digits = 3), 
                round(PseudoR2(fit.h1.decree.nay, "McFadden"), digits = 3), 
                round(PseudoR2(fit.h1.budget.yea, "McFadden"), digits = 3),
                round(PseudoR2(fit.h1.budget.nay, "McFadden"), digits = 3),
                round(PseudoR2(fit.h1.system.yea, "McFadden"), digits = 3),
                round(PseudoR2(fit.h1.system.nay, "McFadden"), digits = 3),
                round(PseudoR2(fit.h1.partis.yea, "McFadden"), digits = 3),
                round(PseudoR2(fit.h1.partis.nay, "McFadden"), digits = 3),
                round(PseudoR2(fit.h1.noretu.yea, "McFadden"), digits = 3),
                round(PseudoR2(fit.h1.noretu.nay, "McFadden"), digits = 3))

# Table 2 in the main paper
stargazer(fit.h1.decree.yea, fit.h1.decree.nay, fit.h1.budget.yea, fit.h1.budget.nay, 
          fit.h1.system.yea, fit.h1.system.nay,
          fit.h1.partis.yea, fit.h1.partis.nay, fit.h1.noretu.yea, fit.h1.noretu.nay,
          se = c(se.h1.decree.yea, se.h1.decree.nay, se.h1.budget.yea, se.h1.budget.nay, 
                 sds, sds2,
                 se.h1.partis.yea, se.h1.partis.nay, se.h1.noretu.yea, se.h1.noretu.nay),
          type="text",
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes        = c("* p<0.1, ** p<0.05, *** p <0.01", 
                           "Clustering-robust standard errors are presented in parentheses.",
                           "For the ordered logistic regression, clustering-robust standard errors are calculated via block-bootstrapping (5000 iterations of resampling)."),
          notes.append = FALSE,
          notes.align = "l",
          dep.var.labels = c("Decree power", "Budget power", "System endorsement", "Partisan president", "No return to parl. system"),
          covariate.labels = c("Treatment: Incumbent loses", "Will vote for the incumbent", "Issue knowledge", "Rightism"),
          column.labels =  c("Initial yea-sayers",  "Initial nay-sayers", "Initial yea-sayers",  "Initial nay-sayers",
                            "Initial yea-sayers",  "Initial nay-sayers","Initial yea-sayers",  "Initial nay-sayers",
                            "Initial yea-sayers",  "Initial nay-sayers"),
          table.layout = "=ldmc#-t-sa=",
          add.lines = list(pseudo.r2s.h1),
          keep.stat = c("n", "ll"),
          digits=3)


# There can be slight differences in the standard errors for Columns 5 and 6 due to the random nature of resampling.
# Log Likelihood  for Columns 5 and 6 are entered manually. 
# For column 5: -605.87 
# For column 6: -378.14 
# (stargazer do not take them to the output)










#### ANALYSES HYPOTHESIS 2: SOCIAL DISTANCE ####
fit.h2.decree.yea.a <- glm(decree_power ~ treatment_erdoganloss * socpol_additive, data=execaggr_yes, family=binomial(link = "logit"))
fit.h2.decree.yea.b <- glm(decree_power ~ treatment_erdoganloss * socpol_additive + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h2.decree.nay.a <- glm(decree_power ~ treatment_erdoganloss * socpol_additive, data=execaggr_no, family=binomial(link = "logit"))
fit.h2.decree.nay.b <- glm(decree_power ~ treatment_erdoganloss * socpol_additive + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

se.h2.decree.yea.a <- list(coef_test(fit.h2.decree.yea.a, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h2.decree.yea.b <- list(coef_test(fit.h2.decree.yea.b, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h2.decree.nay.a <- list(coef_test(fit.h2.decree.nay.a, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.h2.decree.nay.b <- list(coef_test(fit.h2.decree.nay.b, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])

pseudo.r2s.h2 <- c("Psuedo R-sq.", 
                   round(PseudoR2(fit.h2.decree.yea.a, "McFadden"), digits = 3), 
                   round(PseudoR2(fit.h2.decree.yea.b, "McFadden"), digits = 3), 
                   round(PseudoR2(fit.h2.decree.nay.a, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h2.decree.nay.b, "McFadden"), digits = 3))

# Table 4 in the paper
stargazer(fit.h2.decree.yea.a, fit.h2.decree.yea.b,
          fit.h2.decree.nay.a, fit.h2.decree.nay.b,
          se = c(se.h2.decree.yea.a, se.h2.decree.yea.b, 
                 se.h2.decree.nay.a, se.h2.decree.nay.b),
          column.labels =  c("Initial yea-sayers", "Initial yea-sayers", "Initial nay-sayers", "Initial nay-sayers"),
          dep.var.labels = c("Decree power"),
          covariate.labels = c("Treatment: Incumbent loses", "Social distance", "Will vote for the incumbent", "Issue knowledge", "Rightism",
                               "Interaction: Treatment * Social distance" ),
          type="text",
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes        = c("* p<0.1, ** p<0.05, *** p <0.01", 
                           "Clustering-robust standard errors are presented in parentheses."),
          notes.append = FALSE,
          notes.align = "l",
          keep.stat = c("n", "ll"),
          add.lines = list(pseudo.r2s.h2),
          table.layout = "=ldmc#-t-sa=",
          digits=3)

# Rerunning the models for the plots
fit.h2.decree.yea.b <- glm(decree_power ~ tretman * socpol_additive * knowledge + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h2.decree.nay.b <- glm(decree_power ~ tretman * socpol_additive * knowledge + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

# Figure 4
plot_model(fit.h2.decree.yea.b, type = "pred", terms = c("tretman", "socpol_additive[1,4,9]"), 
           title="Initial yea-sayers", axis.title=c("Treatment", "Support for decree power"), axis.lim=c(0.3,1))
plot_model(fit.h2.decree.nay.b, type = "pred", terms = c("tretman", "socpol_additive[1,4,9]"), 
           title="Initial nay-sayers", axis.title=c("Treatment", "Support for decree power"), axis.lim=c(0,0.7))
# The plots are then flipped and processed in a different software. 

#### ANALYSES HYPOTHESIS 3: ECONOMIC ANXIETY ####
fit.h3.decree.yea.a <- glm(decree_power ~ tretman * anxiety_economy, data=execaggr_yes, family=binomial(link = "logit"))
fit.h3.decree.yea.b <- glm(decree_power ~ tretman * anxiety_economy + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h3.decree.nay.a <- glm(decree_power ~ tretman * anxiety_economy, data=execaggr_no, family=binomial(link = "logit"))
fit.h3.decree.nay.b <- glm(decree_power ~ tretman * anxiety_economy + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

se.h3.decree.yea.a <- list(coef_test(fit.h3.decree.yea.a, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h3.decree.yea.b <- list(coef_test(fit.h3.decree.yea.b, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h3.decree.nay.a <- list(coef_test(fit.h3.decree.nay.a, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.h3.decree.nay.b <- list(coef_test(fit.h3.decree.nay.b, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])

pseudo.r2s.h3 <- c("Psuedo R-sq.", 
                   round(PseudoR2(fit.h3.decree.yea.a, "McFadden"), digits = 3), 
                   round(PseudoR2(fit.h3.decree.yea.b, "McFadden"), digits = 3), 
                   round(PseudoR2(fit.h3.decree.nay.a, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h3.decree.nay.b, "McFadden"), digits = 3))

# Table 5 in the paper
stargazer(fit.h3.decree.yea.a, fit.h3.decree.yea.b,
          fit.h3.decree.nay.a, fit.h3.decree.nay.b,
          se = c(se.h3.decree.yea.a, se.h3.decree.yea.b, 
                 se.h3.decree.nay.a, se.h3.decree.nay.b),
          column.labels =  c("Initial yea-sayers", "Initial yea-sayers", "Initial nay-sayers", "Initial nay-sayers"),
          dep.var.labels = c("Decree power"),
          covariate.labels = c("Treatment: Incumbent loses", "Post-incumbent economic anxiety", "Will vote for the incumbent", "Issue knowledge", "Rightism",
                               "Interaction: Treatment * Post-incumbent economic anxiety" ),
          type="text",
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes        = c("* p<0.1, ** p<0.05, *** p <0.01", 
                           "Clustering-robust standard errors are presented in parentheses."),
          notes.append = FALSE,
          notes.align = "l",
          keep.stat = c("n", "ll"),
          add.lines = list(pseudo.r2s.h3),
          table.layout = "=ldmc#-t-sa=",
          digits=3)

# Rerunning the models for the plots
fit.h3.decree.yea.b <- glm(decree_power ~ tretman * anxiety_economy * knowledge + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h3.decree.nay.b <- glm(decree_power ~ tretman * anxiety_economy * knowledge + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

# Figure 5 in the paper
plot_model(fit.h3.decree.yea.b, type = "pred", terms = c("tretman", "anxiety_economy [-2,0,2]"), 
           title="Initial yea-sayers", axis.title=c("Treatment", "Support for decree power"), axis.lim=c(0.3,1))
plot_model(fit.h3.decree.nay.b, type = "pred", terms = c("tretman", "anxiety_economy [-2,0,2]"), 
           title="Initial nay-sayers", axis.title=c("Treatment", "Support for decree power"), axis.lim=c(0,0.7))


#### APPENDIX SECTION 2: EXPLORING WEAKENING TREATMENT EFFECT  ####
# Models
fit.weaken.decree.yea <- glm(decree_power ~ tretman * socpol_additive + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.weaken.decree.nay <- glm(decree_power ~ tretman * socpol_additive  + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))
fit.weaken.budget.yea <- glm(budget_power ~ tretman * socpol_additive  + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.weaken.budget.nay <- glm(budget_power ~ tretman * socpol_additive  + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))
fit.weaken.partis.yea <- glm(executive_president ~ tretman * socpol_additive  + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.weaken.partis.nay <- glm(executive_president ~ tretman * socpol_additive  + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))
fit.weaken.noretu.yea <- glm(noreturn_reduced ~ tretman * socpol_additive  + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.weaken.noretu.nay <- glm(noreturn_reduced ~ tretman * socpol_additive  + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

# Clustering-robust standard errors
se.weaken.decree.yea <- list(coef_test(fit.weaken.decree.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.weaken.decree.nay <- list(coef_test(fit.weaken.decree.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.weaken.budget.yea <- list(coef_test(fit.weaken.budget.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.weaken.budget.nay <- list(coef_test(fit.weaken.budget.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.weaken.partis.yea <- list(coef_test(fit.weaken.partis.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.weaken.partis.nay <- list(coef_test(fit.weaken.partis.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.weaken.noretu.yea <- list(coef_test(fit.weaken.noretu.yea, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.weaken.noretu.nay <- list(coef_test(fit.weaken.noretu.nay, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])


pseudo.r2s.weaken.yea <- c("Psuedo R-sq.", 
                   round(PseudoR2(fit.weaken.decree.yea, "McFadden"), digits = 3), 
                   round(PseudoR2(fit.weaken.budget.yea, "McFadden"), digits = 3),
                   round(PseudoR2(fit.weaken.partis.yea, "McFadden"), digits = 3),
                   round(PseudoR2(fit.weaken.noretu.yea, "McFadden"), digits = 3))

pseudo.r2s.weaken.nay <- c("Psuedo R-sq.", 
                           round(PseudoR2(fit.weaken.decree.nay, "McFadden"), digits = 3), 
                           round(PseudoR2(fit.weaken.budget.nay, "McFadden"), digits = 3),
                           round(PseudoR2(fit.weaken.partis.nay, "McFadden"), digits = 3),
                           round(PseudoR2(fit.weaken.noretu.nay, "McFadden"), digits = 3))

stargazer(fit.weaken.decree.yea, fit.weaken.budget.yea, fit.weaken.partis.yea, fit.weaken.noretu.yea,
          se=c(se.weaken.decree.yea, se.weaken.budget.yea, se.weaken.partis.yea, se.weaken.noretu.yea),
          type="text",
          column.labels =  c("Initial yea-sayers", "Initial yea-sayers", "Initial yea-sayers", "Initial yea-sayers"),
          dep.var.labels = c("Decree power", "Budget power", "Partisan president", "No return to parl. system"),
          covariate.labels = c("Treatment: Incumbent loses", "Social distance", "Will vote for the incumbent", "Issue knowledge", "Rightism",
                               "Interaction: Treatment * Social distance" ),
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes        = c("* p<0.1, ** p<0.05, *** p <0.01", 
                           "Clustering-robust standard errors are presented in parentheses."),
          notes.append = FALSE,
          notes.align = "l",
          keep.stat = c("n", "ll"),
          add.lines = list(pseudo.r2s.weaken.yea),
          table.layout = "=ldmc#-t-sa=",
          digits=3)

stargazer(fit.weaken.decree.nay, fit.weaken.budget.nay, fit.weaken.partis.nay, fit.weaken.noretu.nay,
          se=c(se.weaken.decree.nay, se.weaken.budget.nay, se.weaken.partis.nay, se.weaken.noretu.nay),
          type="text",
          column.labels =  c("Initial nay-sayers", "Initial nay-sayers", "Initial nay-sayers", "Initial nay-sayers"),
          dep.var.labels = c("Decree power", "Budget power", "Partisan president", "No return to parl. system"),
          covariate.labels = c("Treatment: Incumbent loses", "Social distance", "Will vote for the incumbent", "Issue knowledge", "Rightism",
                               "Interaction: Treatment * Social distance" ),
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes        = c("* p<0.1, ** p<0.05, *** p <0.01", 
                           "Clustering-robust standard errors are presented in parentheses."),
          notes.append = FALSE,
          notes.align = "l",
          keep.stat = c("n", "ll"),
          add.lines = list(pseudo.r2s.weaken.nay),
          table.layout = "=ldmc#-t-sa=",
          digits=3)



#### APPENDIX SECTION 3: VARIATION IN ECON. ANX. VARIABLE #### 
prop.table(table(execaggr_yes$anxiety_economy))
prop.table(table(execaggr_no$anxiety_economy))

#### APPENDIX SECTION 4: HYPOTHESIS 2 WITH POLITICAL POLARIZATION ####

fit.h2.decree.yea.a.alt <- glm(decree_power ~ tretman * pol_polar3_alper, data=execaggr_yes, family=binomial(link = "logit"))
fit.h2.decree.yea.b.alt <- glm(decree_power ~ tretman * pol_polar3_alper + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h2.decree.nay.a.alt <- glm(decree_power ~ tretman * pol_polar3_alper, data=execaggr_no, family=binomial(link = "logit"))
fit.h2.decree.nay.b.alt <- glm(decree_power ~ tretman * pol_polar3_alper + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

se.h2.decree.yea.a.alt <- list(coef_test(fit.h2.decree.yea.a.alt, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h2.decree.yea.b.alt <- list(coef_test(fit.h2.decree.yea.b.alt, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h2.decree.nay.a.alt <- list(coef_test(fit.h2.decree.nay.a.alt, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.h2.decree.nay.b.alt <- list(coef_test(fit.h2.decree.nay.b.alt, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])

pseudo.r2s.alt <- c("Psuedo R-sq.", 
                           round(PseudoR2(fit.h2.decree.yea.a.alt, "McFadden"), digits = 3), 
                           round(PseudoR2(fit.h2.decree.yea.b.alt, "McFadden"), digits = 3),
                           round(PseudoR2(fit.h2.decree.nay.a.alt, "McFadden"), digits = 3),
                           round(PseudoR2(fit.h2.decree.nay.b.alt, "McFadden"), digits = 3))

stargazer(fit.h2.decree.yea.a.alt, fit.h2.decree.yea.b.alt,
          fit.h2.decree.nay.a.alt, fit.h2.decree.nay.b.alt,
          se = c(se.h2.decree.yea.a.alt, se.h2.decree.yea.b.alt, 
                 se.h2.decree.nay.a.alt, se.h2.decree.nay.b.alt),
          type="text",
          column.labels =  c("Initial yea-sayers", "Initial yea-sayers", "Initial nay-sayers", "Initial nay-sayers"),
          dep.var.labels = c("Decree power"),
          covariate.labels = c("Treatment: Incumbent loses", "Partisan distance", "Will vote for the incumbent", "Issue knowledge", "Rightism",
                               "Interaction: Treatment * Partisan distance" ),
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes        = c("* p<0.1, ** p<0.05, *** p <0.01", 
                           "Clustering-robust standard errors are presented in parentheses."),
          notes.append = FALSE,
          notes.align = "l",
          keep.stat = c("n", "ll"),
          add.lines = list(pseudo.r2s.alt),
          table.layout = "=ldmc#-t-sa=",
          digits=3)






#### APPENDIX SECTION 5.1: DEMOGRAPHIC CORRELATES OF MODERATORS ####
fit.a <- lm(data=execaggr_yes, socpol_additive ~ age + male + edu_level + income_household + sunni + kurd)
fit.b <- lm(data=execaggr_yes, anxiety_economy ~ age + male + edu_level + income_household + sunni + kurd)
fit.c <- lm(data=execaggr_no, socpol_additive ~ age + male + edu_level + income_household + sunni + kurd)
fit.d <- lm(data=execaggr_no, anxiety_economy ~ age + male + edu_level + income_household + sunni + kurd)


stargazer(fit.a, fit.b, fit.c, fit.d,
          type="text",
          dep.var.labels = c("Social distance", "Post-incumbent economic anxiety", "Social distance", "Post-incumbent economic anxiety"),
          covariate.labels = c("Age", "Male", "Education level", "Household income",
                               "Sunni", "Kurdish"),
          column.labels = c("Initial yea-sayers", "Initial yea-sayers", "Initial nay-sayers","Initial nay-sayers"),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Normal standard errors are presented in parentheses."))






#### APPENDIX SECTION 5.2: BOTH MODERATORS IN THE MODEL ####

# What is the correlation between economic anxiety and social distance?
cor(execaggr_yes$anxiety_economy, execaggr_yes$socpol_additive, use="complete.obs", method="pearson")
cor(execaggr_no$anxiety_economy, execaggr_no$socpol_additive, use="complete.obs", method="pearson")

fit.h2h3.decree.yea.a <- glm(decree_power ~ tretman * socpol_additive + tretman * anxiety_economy, data=execaggr_yes, family=binomial(link = "logit"))
fit.h2h3.decree.yea.b <- glm(decree_power ~ tretman * socpol_additive + tretman * anxiety_economy + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h2h3.decree.nay.a <- glm(decree_power ~ tretman * socpol_additive + tretman * anxiety_economy, data=execaggr_no, family=binomial(link = "logit"))
fit.h2h3.decree.nay.b <- glm(decree_power ~ tretman * socpol_additive + tretman * anxiety_economy + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

se.h2h3.decree.yea.a <- list(coef_test(fit.h2h3.decree.yea.a, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h2h3.decree.yea.b <- list(coef_test(fit.h2h3.decree.yea.b, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h2h3.decree.nay.a <- list(coef_test(fit.h2h3.decree.nay.a, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])
se.h2h3.decree.nay.b <- list(coef_test(fit.h2h3.decree.nay.b, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])

pseudo.r2s.h2h3 <- c("Psuedo R-sq.", 
                    round(PseudoR2(fit.h2h3.decree.yea.a, "McFadden"), digits = 3), 
                    round(PseudoR2(fit.h2h3.decree.yea.b, "McFadden"), digits = 3),
                    round(PseudoR2(fit.h2h3.decree.nay.a, "McFadden"), digits = 3),
                    round(PseudoR2(fit.h2h3.decree.nay.b, "McFadden"), digits = 3))

stargazer(fit.h2h3.decree.yea.a, fit.h2h3.decree.yea.b,
          fit.h2h3.decree.nay.a, fit.h2h3.decree.nay.b,
          se = c(se.h2h3.decree.yea.a, se.h2h3.decree.yea.b, 
                 se.h2h3.decree.nay.a, se.h2h3.decree.nay.b),
          column.labels =  c("Initial yea-sayers", "Initial yea-sayers", "Initial nay-sayers", "Initial nay-sayers"),
          dep.var.labels = c("Decree power"),
          covariate.labels = c("Treatment: Incumbent loses", "Social distance", "Post-incumbent economic anxiety", 
                               "Will vote for the incumbent", "Issue knowledge", "Rightism",
                               "Interaction: Treatment * Social distance", "Interaction: Treatment * Post-incumbent economic anxiety" ),
          type="text",
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes        = "*p<0.1; **p<0.05; ***p<0.01. Clustering-robust standard errors are presented parentheses.", 
          notes.append = FALSE,
          notes.align = "l",
          keep.stat = c("n", "ll"),
          add.lines = list(pseudo.r2s.h2h3),
          table.layout = "=ldmc#-t-sa=",
          digits=3)

#### APPENDIX SECTION 8: ADDITIONAL ANALYSES WITH ALTERNATIVE DATASETS ####
## Dividing the datasets for the additional analyses
execaggr_proerdogan <- execaggr[(execaggr$vote_firstround == 1 & is.na(execaggr$vote_firstround) == 0),]
execaggr_antierdogan <- execaggr[((execaggr$vote_firstround == 2 | execaggr$vote_firstround == 3 | 
                                 execaggr$vote_firstround == 4 | execaggr$vote_firstround == 5) & is.na(execaggr$vote_firstround) == 0),]

execaggr_procumhur2015 <- execaggr[(execaggr$past_vote == 1 | execaggr$past_vote == 4),]
execaggr_proopp2015 <- execaggr[(execaggr$past_vote == 2 | execaggr$past_vote == 3),]


# Analyses based on who they voted in 2015 elections
execaggr_procumhur2015$tretman <- ifelse(execaggr_procumhur2015$treatment_erdoganloss==1, "2 Incumbent Loses", "1 Incumbent Wins")
execaggr_proopp2015$tretman <- ifelse(execaggr_proopp2015$treatment_erdoganloss==1, "2 Incumbent Loses", "1 Incumbent Wins")

# Models
fit.h1.decree.yea.alt1 <- glm(decree_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_procumhur2015, family=binomial(link = "logit"))
fit.h1.decree.nay.alt1 <- glm(decree_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_proopp2015, family=binomial(link = "logit"))
fit.h1.budget.yea.alt1 <- glm(budget_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_procumhur2015, family=binomial(link = "logit"))
fit.h1.budget.nay.alt1 <- glm(budget_power ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_proopp2015, family=binomial(link = "logit"))
fit.h1.system.yea.alt1 <- polr(factor(system_endorsement) ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_procumhur2015, Hess=T)
fit.h1.system.nay.alt1 <- polr(factor(system_endorsement) ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_proopp2015, Hess=T)
fit.h1.partis.yea.alt1 <- glm(executive_president ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_procumhur2015, family=binomial(link = "logit"))
fit.h1.partis.nay.alt1 <- glm(executive_president ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_proopp2015, family=binomial(link = "logit"))
fit.h1.noretu.yea.alt1 <- glm(noreturn_reduced ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_procumhur2015, family=binomial(link = "logit"))
fit.h1.noretu.nay.alt1 <- glm(noreturn_reduced ~ tretman + rte_firstround + knowledge + left_right, data=execaggr_proopp2015, family=binomial(link = "logit"))

# Clustering-robust standard errors
se.h1.decree.yea.alt1 <- list(coef_test(fit.h1.decree.yea.alt1, vcov = "CR2", cluster = execaggr_procumhur2015$ilcemah, test = "Satterthwaite")[,2])
se.h1.decree.nay.alt1 <- list(coef_test(fit.h1.decree.nay.alt1, vcov = "CR2", cluster = execaggr_proopp2015$ilcemah, test = "Satterthwaite")[,2])
se.h1.budget.yea.alt1 <- list(coef_test(fit.h1.budget.yea.alt1, vcov = "CR2", cluster = execaggr_procumhur2015$ilcemah, test = "Satterthwaite")[,2])
se.h1.budget.nay.alt1 <- list(coef_test(fit.h1.budget.nay.alt1, vcov = "CR2", cluster = execaggr_proopp2015$ilcemah, test = "Satterthwaite")[,2])
se.h1.system.yea.alt1 <- list(summary(fit.h1.system.yea.alt1)$coefficients[,2])
se.h1.system.nay.alt1 <- list(summary(fit.h1.system.nay.alt1)$coefficients[,2])
se.h1.partis.yea.alt1 <- list(coef_test(fit.h1.partis.yea.alt1, vcov = "CR2", cluster = execaggr_procumhur2015$ilcemah, test = "Satterthwaite")[,2])
se.h1.partis.nay.alt1 <- list(coef_test(fit.h1.partis.nay.alt1, vcov = "CR2", cluster = execaggr_proopp2015$ilcemah, test = "Satterthwaite")[,2])
se.h1.noretu.yea.alt1 <- list(coef_test(fit.h1.noretu.yea.alt1, vcov = "CR2", cluster = execaggr_procumhur2015$ilcemah, test = "Satterthwaite")[,2])
se.h1.noretu.nay.alt1 <- list(coef_test(fit.h1.noretu.nay.alt1, vcov = "CR2", cluster = execaggr_proopp2015$ilcemah, test = "Satterthwaite")[,2])

pseudo.r2s.h1.alt1 <- c("Psuedo R-sq.", 
                   round(PseudoR2(fit.h1.decree.yea.alt1, "McFadden"), digits = 3), 
                   round(PseudoR2(fit.h1.decree.nay.alt1, "McFadden"), digits = 3), 
                   round(PseudoR2(fit.h1.budget.yea.alt1, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h1.budget.nay.alt1, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h1.system.yea.alt1, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h1.system.nay.alt1, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h1.partis.yea.alt1, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h1.partis.nay.alt1, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h1.noretu.yea.alt1, "McFadden"), digits = 3),
                   round(PseudoR2(fit.h1.noretu.nay.alt1, "McFadden"), digits = 3))


stargazer(fit.h1.decree.yea.alt1, fit.h1.decree.nay.alt1, fit.h1.budget.yea.alt1, fit.h1.budget.nay.alt1, fit.h1.system.yea.alt1, fit.h1.system.nay.alt1,
          fit.h1.partis.yea.alt1, fit.h1.partis.nay.alt1, fit.h1.noretu.yea.alt1, fit.h1.noretu.nay.alt1,
          se = c(se.h1.decree.yea.alt1, se.h1.decree.nay.alt1, se.h1.budget.yea.alt1, se.h1.budget.nay.alt1, se.h1.system.yea.alt1, se.h1.system.nay.alt1,
                 se.h1.partis.yea.alt1, se.h1.partis.nay.alt1, se.h1.noretu.yea.alt1, se.h1.noretu.nay.alt1),
          dep.var.labels = c("Decree power", "Budget power", "System endorsement", "Partisan president", "No return to parl. system"),
          covariate.labels = c("Treatment: Incumbent loses", "Will vote for incumbent", "Issue knowledge", "Rightism"),
          column.labels = c("Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp."),
          type="text",
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the neighborhood level) are presented in brackets, except for models for the dependent variable 'System endorsement,' in which normal standard errors are presented."),
          keep.stat = c("n", "ll"),
          add.lines = list(pseudo.r2s.h1.alt1),
          table.layout = "=ldmc#-t-sa=",
          star.cutoffs = c(0.1, 0.05, 0.01), 
          digits=3)

stats::logLik(fit.h1.system.yea.alt1)
stats::logLik(fit.h1.system.nay.alt1)


# Analyses based on who they intend to vote in 2018 elections
execaggr_proerdogan$tretman <- ifelse(execaggr_proerdogan$treatment_erdoganloss==1, "2 Incumbent Loses", "1 Incumbent Wins")
execaggr_antierdogan$tretman <- ifelse(execaggr_antierdogan$treatment_erdoganloss==1, "2 Incumbent Loses", "1 Incumbent Wins")


# Models
fit.h1.decree.yea.alt2 <- glm(decree_power ~ tretman + knowledge + left_right, data=execaggr_proerdogan, family=binomial(link = "logit"))
fit.h1.decree.nay.alt2 <- glm(decree_power ~ tretman + knowledge + left_right, data=execaggr_antierdogan, family=binomial(link = "logit"))
fit.h1.budget.yea.alt2 <- glm(budget_power ~ tretman + knowledge + left_right, data=execaggr_proerdogan, family=binomial(link = "logit"))
fit.h1.budget.nay.alt2 <- glm(budget_power ~ tretman + knowledge + left_right, data=execaggr_antierdogan, family=binomial(link = "logit"))
fit.h1.system.yea.alt2 <- polr(factor(system_endorsement) ~ tretman + knowledge + left_right, data=execaggr_proerdogan, Hess=T)
fit.h1.system.nay.alt2 <- polr(factor(system_endorsement) ~ tretman + knowledge + left_right, data=execaggr_antierdogan, Hess=T)
fit.h1.partis.yea.alt2 <- glm(executive_president ~ tretman + knowledge + left_right, data=execaggr_proerdogan, family=binomial(link = "logit"))
fit.h1.partis.nay.alt2 <- glm(executive_president ~ tretman + knowledge + left_right, data=execaggr_antierdogan, family=binomial(link = "logit"))
fit.h1.noretu.yea.alt2 <- glm(noreturn_reduced ~ tretman + knowledge + left_right, data=execaggr_proerdogan, family=binomial(link = "logit"))
fit.h1.noretu.nay.alt2 <- glm(noreturn_reduced ~ tretman + knowledge + left_right, data=execaggr_antierdogan, family=binomial(link = "logit"))

# Clustering-robust standard errors
se.h1.decree.yea.alt2 <- list(coef_test(fit.h1.decree.yea.alt2, vcov = "CR2", cluster = execaggr_proerdogan$ilcemah, test = "Satterthwaite")[,2])
se.h1.decree.nay.alt2 <- list(coef_test(fit.h1.decree.nay.alt2, vcov = "CR2", cluster = execaggr_antierdogan$ilcemah, test = "Satterthwaite")[,2])
se.h1.budget.yea.alt2 <- list(coef_test(fit.h1.budget.yea.alt2, vcov = "CR2", cluster = execaggr_proerdogan$ilcemah, test = "Satterthwaite")[,2])
se.h1.budget.nay.alt2 <- list(coef_test(fit.h1.budget.nay.alt2, vcov = "CR2", cluster = execaggr_antierdogan$ilcemah, test = "Satterthwaite")[,2])
se.h1.system.yea.alt2 <- list(summary(fit.h1.system.yea.alt2)$coefficients[,2])
se.h1.system.nay.alt2 <- list(summary(fit.h1.system.nay.alt2)$coefficients[,2])
se.h1.partis.yea.alt2 <- list(coef_test(fit.h1.partis.yea.alt2, vcov = "CR2", cluster = execaggr_proerdogan$ilcemah, test = "Satterthwaite")[,2])
se.h1.partis.nay.alt2 <- list(coef_test(fit.h1.partis.nay.alt2, vcov = "CR2", cluster = execaggr_antierdogan$ilcemah, test = "Satterthwaite")[,2])
se.h1.noretu.yea.alt2 <- list(coef_test(fit.h1.noretu.yea.alt2, vcov = "CR2", cluster = execaggr_proerdogan$ilcemah, test = "Satterthwaite")[,2])
se.h1.noretu.nay.alt2 <- list(coef_test(fit.h1.noretu.nay.alt2, vcov = "CR2", cluster = execaggr_antierdogan$ilcemah, test = "Satterthwaite")[,2])

pseudo.r2s.h1.alt2 <- c("Psuedo R-sq.", 
                        round(PseudoR2(fit.h1.decree.yea.alt2, "McFadden"), digits = 3), 
                        round(PseudoR2(fit.h1.decree.nay.alt2, "McFadden"), digits = 3), 
                        round(PseudoR2(fit.h1.budget.yea.alt2, "McFadden"), digits = 3),
                        round(PseudoR2(fit.h1.budget.nay.alt2, "McFadden"), digits = 3),
                        round(PseudoR2(fit.h1.system.yea.alt2, "McFadden"), digits = 3),
                        round(PseudoR2(fit.h1.system.nay.alt2, "McFadden"), digits = 3),
                        round(PseudoR2(fit.h1.partis.yea.alt2, "McFadden"), digits = 3),
                        round(PseudoR2(fit.h1.partis.nay.alt2, "McFadden"), digits = 3),
                        round(PseudoR2(fit.h1.noretu.yea.alt2, "McFadden"), digits = 3),
                        round(PseudoR2(fit.h1.noretu.nay.alt2, "McFadden"), digits = 3))

stargazer(fit.h1.decree.yea.alt2, fit.h1.decree.nay.alt2, fit.h1.budget.yea.alt2, fit.h1.budget.nay.alt2, fit.h1.system.yea.alt2, fit.h1.system.nay.alt2,
          fit.h1.partis.yea.alt2, fit.h1.partis.nay.alt2, fit.h1.noretu.yea.alt2, fit.h1.noretu.nay.alt2,
          se = c(se.h1.decree.yea.alt2, se.h1.decree.nay.alt2, se.h1.budget.yea.alt2, se.h1.budget.nay.alt2, se.h1.system.yea.alt2, se.h1.system.nay.alt2,
                 se.h1.partis.yea.alt2, se.h1.partis.nay.alt2, se.h1.noretu.yea.alt2, se.h1.noretu.nay.alt2),
          dep.var.labels = c("Decree power", "Budget power", "System endorsement", "Partisan president", "No return to parl. system"),
          covariate.labels = c("Treatment: Incumbent loses",  "Issue knowledge", "Rightism"),
          column.labels = c("Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp.", "Pro-gov.", "Pro-opp."),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the neighborhood level) are presented in brackets, except for models for the dependent variable 'System endorsement,' in which normal standard errors are presented."),
          type="text",
          star.cutoffs = c(0.1, 0.05, 0.01), 
          add.lines = list(pseudo.r2s.h1.alt2),
          table.layout = "=ldmc#-t-sa=",
          keep.stat = c("n", "ll"),
          digits=3)

stats::logLik(fit.h1.system.yea.alt2)
stats::logLik(fit.h1.system.nay.alt2)

#### APPENDIX SECTION 9: PARTISAN PRESIDENT AND RIGHT-WING IDEOLOGY ####
fit.h1.exec.nay.right <- glm(executive_president ~ treatment_erdoganloss * left_right + rte_firstround + knowledge, data=execaggr_no, family=binomial(link = "logit"))
se.h1.exec.nay.right <- list(coef_test(fit.h1.exec.nay.right, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])

pseudo.r2s.right <- c("Psuedo R-sq.", 
                        round(PseudoR2(fit.h1.exec.nay.right, "McFadden"), digits = 3))

stargazer(fit.h1.exec.nay.right, 
          se = c(se.h1.exec.nay.right),
          type="text",
          covariate.labels = c("Treatment: Incumbent loses", "Rightism",
                               "Will vote for the incumbent", "Issue knowledge", "Rightism", "Treatmment * Rightism"),
          dep.var.labels = "Partisan president",
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the neighborhood level) are presented in parentheses."),
          add.lines = list(pseudo.r2s.right),
          table.layout = "=ldmc#-t-sa=",
          keep.stat = c("n", "ll"),
          digits=3)



#### APPENDIX SECTION 10: PRIOR EXPECTATION #### 

# Models
fit.h1.decree.yea.prior <- glm(decree_power ~ tretman *  prior_erdoganwin + rte_firstround + knowledge + left_right, data=execaggr_yes, family=binomial(link = "logit"))
fit.h1.decree.nay.prior <- glm(decree_power ~ tretman *  prior_erdoganwin + rte_firstround + knowledge + left_right, data=execaggr_no, family=binomial(link = "logit"))

# Clustering-robust standard errors
se.h1.decree.yea.prior <- list(coef_test(fit.h1.decree.yea.prior, vcov = "CR2", cluster = execaggr_yes$ilcemah, test = "Satterthwaite")[,2])
se.h1.decree.nay.prior <- list(coef_test(fit.h1.decree.nay.prior, vcov = "CR2", cluster = execaggr_no$ilcemah, test = "Satterthwaite")[,2])


pseudo.r2s.h1.prior <- c("Psuedo R-sq.", 
                        round(PseudoR2(fit.h1.decree.yea.prior, "McFadden"), digits = 3), 
                        round(PseudoR2(fit.h1.decree.nay.prior, "McFadden"), digits = 3))

stargazer(fit.h1.decree.yea.prior, fit.h1.decree.nay.prior, 
          se = c(se.h1.decree.yea.prior, se.h1.decree.nay.prior),
          type="text",
          covariate.labels = c("Treatment: Incumbent loses", "Prior expectation: Incumbent win", 
                               "Will vote for the incumbent", "Issue knowledge", "Rightism", "Treatmment * Prior expectation"),
          dep.var.labels = "Support for decree power",
          column.labels = c("Initial yea-sayers", "Initial nay-sayers"),
          star.cutoffs = c(0.1, 0.05, 0.01), 
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the neighborhood level) are presented in parentheses."),
          add.lines = list(pseudo.r2s.h1.prior),
          table.layout = "=ldmc#-t-sa=",
          keep.stat = c("n", "ll"),
          digits=3)

